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Abstract. The method of recurrence plots is extended to the 
cross recurrence plots (CRP), which among others enables 
the study of synchronization or time differences in two time 
series. This is emphasized in a distorted main diagonal in 
the cross recurrence plot, the Une of synchronization (LOS). 
A non-parametrical fit of this LOS can be used to rescale 
the time axis of the two data series (whereby one of it is 
e. g. compressed or stretched) so that they are synchronized. 
An application of this method to geophysical sediment core 
data illustrates its suitability for real data. The rock magnetic 
data of two different sediment cores from the Makarov Basin 
can be adjusted to each other by using this method, so that 
they are comparable. 



1 Introduction 

The adjustment of data sets with various time scales occurs 
in many occasions, e. g. data preparation of tree rings or geo- 
physical profiles. In geology, often a large set of geophysi- 
cal data series is gained at various locations (e. g. sediment 
cores). That is why these data series have a different length 
and time scale. Before any time series analysis can be started, 
the data series have to be synchronized to the same time 
scale. Usually, this is done visually by comparing and cor- 
relating each maximum and minimum in both data sets by 
hand (wiggle matching), which includes the human factor of 
subjectiveness and is a lengthy process. An automatical and 
objective method for verification should be very welcome. 

In the last decades some techniques for this kind of corre- 
lation and adjustment were suggested. They span graphical 
methods (Prell et al., 1986), inverse algorithms, e.g. using 
Fourier series (Martinson et al., 1982) and algorithms based 
on similarity of data, e. g. sequence slotting (Thompson and 
Clark, 1989). 

However, we focus on a method based on nonlinear time 
series analysis. During our investigations of the method of 
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cross recurrence plots (CRP), we have found an interesting 
feature of it. Besides the possibiUty of the application of the 
recurrence quantification analysis (RQA) of Webber and Zbi- 
lut on CRPs (1998), there is a more fundamental relation be- 
tween the structures in the CRP and the considered systems. 
Finally, this feature can be used for the task of the synchro- 
nization of data sets. Although the first steps of this method 
are similar to the sequence slotting method, their roots are 
different. 

First we give an introduction in CRPs. Then we explain 
the relationship between the structures in the CRP and the 
systems and illustrate this with a simple model. Finally we 
apply the CRP on geophysical data in order to synchronize 
various profiles and to show their practical availability. Since 
we focus on the synchronization feature of the CRP, we will 
not give a comparison between the different alignment meth- 
ods. 



2 The Recurrence Plot 

Recurrence plots (RP) were firstly introduced by Eckmann 
et al. (1987) in order to visualize time dependent behaviour 
of orbits Xj in the phase space. A RP represents the recur- 
rence of the phase space trajectory to a state. The recurrence 
of states is a fundamental property of deterministic dynami- 
cal systems (Argyris et al., 1994; Casdagli, 1997; Kantz and 
Schreiber, 1997). The main step of the visualization is the 
calculation of the N x A'^-matrix 

Tiij = @{e-\\xi-Xj\\), i,j = l...N, (1) 

where e is a predefined cut-off distance, || • || is the norm 
(e.g. the Euchdean norm) and 6(x) is the Heaviside func- 
tion. The values one and zero in this matrix can be simply 
visualized by the colours black and white. Depending on the 
kind of the apphcation, e can be a fixed value or it can be 
changed for each i in such a way that in the ball with the ra- 
dius s a predefined amount of neighbours occurs. The latter 
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will provide a constant density of recurrence points in each 
column of the RP. 

The recurrence plot exhibits characteristic patterns for typ- 
ical dynamical behaviour (Eckmann et al., 1987; Webber Jr 
and Zbilut, 1994): A collection of single recurrence points, 
homogeneously and irregularly distributed over the whole 
plot, reveals a mainly stochastic process. Longer, parallel 
diagonals formed by recurrence points and with the same 
distance between the diagonals are caused by periodic pro- 
cesses. A paling of the RP away from the main diagonal to 
the corners reveals a drift in the amplitude of the system. 
Vertical and horizontal white bands in the RP result from 
states which occur rarely or represent extreme states. Ex- 
tended horizontal and vertical black lines or areas occur if a 
state does not change for some time, e. g. laminar states. All 
these structures were formed by using the property of recur- 
rence of states. It should be pointed out that the states are 
only the "same" and recur in the sense of the vicinity, which 
is determined by the distance e. RPs and their quantitative 
analysis (RQA) became more well known in the last decade 
(e. g. Casdagli, 1997). Their applications to a wide field of 
miscellaneous research show their suitability in the analysis 
of short and non-stationary data. 

3 The Cross Recurrence Plot 

Analogous to Zbilut et al. (1998), we have expanded the 
method of recurrence plots (RP) to the method of cross re- 
currence plots. In contrast to the conventional RP, two time 
series are simultaneously embedded in the same phase space. 
The test for closeness of each point of the first trajectory Xi 
(i = 1 . . . N) with each point of the second trajectory yj 
{] = !... M) results in a x M array 

CR„ =e(£-||x, -^,11). (2) 

Its visualization is called the cross recurrence plot (CRP). 
The definition of the closeness between both trajectories can 
be varied as described above. Varying e may be useful to 
handle systems with different amplitudes. 

The CRP compares the considered systems and allows us 
to benchmark the similarity of states. In this paper we fo- 
cus on the bowed ''main diagonal" in the CRP, because it is 
related to the frequencies and phases of the systems consid- 
ered. 

4 The Line of Synchronization in the CRP 

Regarding the conventional RP, Eq. 1, one always finds a 
main diagonal in the plot, because of the identity of the (i, i)- 
states. The RP can be considered as a special case of the 
CRP, Eq. 2, which usually does not have a main diagonal as 
the (i, z)-states are not identical. 

In data analysis one is often faced with time series that 
are measured on varying time scales. These could be sets 
from borehole or core data in geophysics or tree rings in 
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Fig. 1. Cross recurrence plots of sine functions /(i) = sin((/3t) and g(t) = 
sin(ipt + asin(i/>t)), whereat a = for the black CRP, a = 0.5 for the 
green CRP and a = 1 for the red CRP. The variation in the time domain 
leads to a deforming of the synchronization line. 

dendrochronology. Sediment cores might have undergone a 
number of coring disturbances such as compression or stretch- 
ing. Moreover, cores from different sites with differing sed- 
imentation rates would have different temporal resolutions. 
All these factors require a method of synchronizing. 

A CRP of the two corresponding time series will not con- 
tain a main diagonal. But if the sets of data are similar 
e. g. only rescaled, a more or less continuous line in the CRP 
that is like a distorted main diagonal can occur. This line 
contains information on the rescaling. We give an illustrative 
example. A CRP of a sine function with itself (i. e. this is the 
RP) contains a main diagonal (black CRP in Fig. 1). Hence, 
the CRPs in the Fig. 1 are computed with embeddings of di- 
mension one, further diagonal lines from the upper left to 
the lower right occur. These lines typify the similarity of the 
phase space trajectories in positive and negative time direc- 
tion. 

Now we rescale the time axis of the second sine function 
in the following way 

sin(iy9t) — > sin((y9t + a sin('0i)) (3) 

We will henceforth use the notion rescaling only in the 
mention of the rescaling of the time scale. The rescaling of 
the second sine function with different parameters Lp results 
in a deformation of the main diagonal (green and red CRP 
in Fig. 1). The distorted line contains the information on the 
rescaling, which we will need in order to re-synchronize the 
two time series. Therefore, we call this distorted diagonal 
line of synchronization (LOS). 

In the following, we present a toy function to explain the 
procedure. If we consider a one dimensional case without 
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embedding, the CRP is computed with 
CR{h,t2) = &{e-\\f{h)-git2)\\). 



(4) 



If we set e = to simplify the condition, Eq. (4) delivers 
a recurrence point if 



fih) = g{t2). 



(5) 



In general, this is an impUcit condition that links the variable 
ti to t2. Considering the physical examples of above, it can 
be assumed that the time series are essentially the same - that 
means that / = 51 - up to a rescaling function of time. So we 
can state 



(6) 



If the functions /(•) and g{-) are not identical our method is 
in general not capable of deciding if the difference in the time 
series is due to different dynamics (/(•) ^ g{-)) or if it is due 
to simple rescaling. So the assumption that the dynamics is 
alike up to a rescaling in time is essential. Even though for 
some cases where f ^ git can be applied in the same way. If 
we consider the functions /(•) = a- f{-) + b and g{-) = §{■), 
whereby /(•) ^ g{-) are the observations and /(•) = g{-) 
are the states, normalization with respect to the mean and the 
standard deviation allows us to use our method. 



/(•) 

K-) 



a-f{-) + b~ 



/(•) = 



/(•) - (/(•)) 
^ (/(•)) 



(7) 
(8) 



With g{-) = /(•) the functions /(•) and g{-) are the same af- 
ter the normalization. Then our method can be applied with- 
out any further modification. 

In some special cases Eq. (6) can be resolved with respect 
to ti. Such a case is a system of two sine functions with 
different frequencies 



m 

9{t) 



sm{ip ■t + a) 
sm{tp -t + p) 



Using Eq. (5) and Eq. (6) we find 

sin {(fi ti + a) — sin {il)t2+ (3) = 

and one explicit solution of this equation is 



i2 = 0(ii)= ( 1^1+7 



(9) 
(10) 

(11) 

(12) 



with 7 = . In this special case the slope of the main 
line in a cross recurrence plot represents the frequency ratio 
and the distance between the axes origin and the intersec- 
tion of the line of synchronization with the ordinate reveals 
the phase difference. The function t2 = 4>{ti) (Eq. 6) is a 
transfer or rescaling function which allows us to rescale the 
second system to the first system. If the rescaling function is 
not hnear the LOS will also be curved. 

For the application one has to determine the LOS - usually 
non-parametrically - and then rescale one of the time series. 



In the Appendix we describe a simple algorithm for esti- 
mating the LOS. Its determination will be better for higher 
embeddings because the vertical and cross-diagonal struc- 
tures will vanish. We do not conceal that the embedding of 
the time series is connected with difficulties. The Takens 
Embedding Theorem holds for closed, deterministic systems 
without noise, only. If noise is present one needs its real- 
ization to find a reasonable embedding. For stochastic time 
series it does not make sense to consider a phase space and 
so embedding is in general not justified here either (Romano, 
to be published; Takens, 1981). 

The choice of a special embedding lag could be correct 
for one section of the data, but incorrect for another (for an 
example see below). This can be the case if the data is non- 
stationary. Furthermore, the choice of method of computing 
the CRP and the threshold e will influence the quality of the 
estimated LOS. 

The next sections will be dedicated to application. 

5 Application to a Simple Example 

At first, we consider two sine functions f{t) = s'm{(pt) and 
g{t) = sm{ipt^), where the time scale of the second sine 
differs from the first by a quadratic term and the frequency 
tjj = 0.01 (fi. Sediment parameters are related to such kind 
of functions because gravity pressure increases nonlinearly 
with the depth. It can be assumed that both data series come 
from the same process and were subjected to different de- 
posital compressions (e. g. a squared or exponential increas- 
ing of the compression). Their CRP contains a bowed LOS 
(Fig. 2). We have used the embedding parameters dimension 
TO = 2, delay r = 7r/2 and a varying threshold s, so that the 
CRP contains a constant recurrence density of 20 %. Assum- 
ing that the time scale of g is not the correct scale, we denote 
that scale by t'. In order to determine the non-parametrical 
LOS, we have implemented the algorithm described in the 
Appendix. Although this algorithm is still not mature, we ob- 
tained rehable results (Fig. 3). The resulting rescaling func- 
tion has the expected squared shape t = (j){t') = 0.01 1'^ 
(red curve in Fig. 3). Substituting the time scale t' in the sec- 
ond data series g(t') by this rescaling function t = 4>{t'),we 
get a set of synchronized data f{t) and g{t) with the non- 
parametric rescaling function t = 4>{t') (Fig. 4). The syn- 
chronized data series are approximately the same. The cause 
of some differences is the meandering of the LOS, which it- 
self is caused by partial weak embedding. Nevertheless, this 
can be avoided by more complex algorithm for estimating the 
LOS. 



6 Application to Real Data 

In order to continue the illustration of the working of our 
method we have applied it to real data from geology. 

In the following we compare the method of cross recur- 
rence plot matching with the conventional method of visual 
wiggle matching (interactive adjustment). Geophysical data 
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Fig. 2. Cross recurrence plots of two sine functions f{t) = sin{ipt) and 
g{t) = sin{rpt^)) wliicli is tlie base for tlie determination of the rescaling 
function between both data series. The embedding parameters were dimen- 
sion m = 2, delay t = n/2 and a varying threshold e, in such a way that 
the CRP contains a constant recuiTence density of 20 %. 



Reseating Function 




Fig. 3. The rescaling function (black) determined from the CRP in Fig. 2. 
It has the expected parabolic shape of the squared coherence in the time 
domain. In red the square function. 




Fig. 4. Reference data series (upper panel) and rescaled data series before 
(red) and after (black) rescaling by using the rescaling function of Fig. 3 
(lower panel). 



of two sediment cores from the Makarov Basin, central Arc- 
tic Ocean, PS 2178-3 and PS 2180-2, were analysed. The 
task should be to adjust the data of the PS 2178-3 data (data 
length N = 436) to the scale of the PS 2180-2 (data length 
N = 251) in order to get a depth-depth-function which al- 
lows to synchronize both data sets (Fig. 5). 

We have constructed the phase space with the normalized 
six parameters low field magnetic susceptibility (klf), an- 
hysteretic remanent magnetization (ARM), ratio of anhys- 
teretic susceptibility to klf (karm / ^^lf), relative palaeoin- 
tensity ( PJA), median destructive field of ARM {MDFarm) 
and inclination (INC). A comprehensive discussion of the 
data is given in Nowaczyk et al. (2001). The embedding was 
combined with the time-delayed method according to (Tak- 
ens, 1981) in order to increase further the dimension of the 
phase-space with the following rule: If we have n parameters 
Gi, the embedding with dimension m and delay r will result 
in a (to • n)-dimensional phase space: 

x{t) = (ai(t), . . . ,a„(t), 

ai{t + T),...,an{t + T), 

ai{t + 2T),...,a^{t + 2T),... 

ai{t + (m - 1)t), . . . , Unit + (to - 1)t) (13) 

For our investigation we have used a dimension m = 3 
and a delay t = 1, which finally led to a phase space of 
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Fig. 5. ARM data of the boreholes PS 2178-3 GPC and PS 2180-2 GPC in 
the Central Arctic Ocean before adjustment. 
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Depth in Core PS 2178-3 [cm] 
200 400 600 800 1000 1200 1400 




Fig. 6. Cross recuiTence plot based on six normalized sediment parameters 
and an additional embedding dimension of m = 3 (r = 1, e = 0.05). 



dimension 18 (3 x 6). The recurrence criterion was e = 5 % 
nearest neighbours. 

The resulting CRP shows a clear LOS and some clustering 
of black patches (Fig. 6). The latter occurs due to the plateaus 
in the data. The next step is to fit a non-parametric function 
(the depth-depth-curve) to the LOS in the CRP (red curve in 
Fig. 6). With this function we are able to adjust the data of 
the PS 2178-3 core to the scale of PS 2180-2 (Fig. 8). 

The determination of the depth-depth-function with the 
conventional method of visual wiggle matching is based on 
the interactive and parallel searching for the same structures 
in the different parameters of both data sets. If the adjust- 
ment does not work in a section of the one parameter, one 
can use another parameter for this section, which allows the 
multivariate adjustment of the data sets. The recognition of 
the same structures in the data sets requires a degree of expe- 
rience. However, human eyes are usually better in the visual 
assessment of complex structures than a computational algo- 
rithm. 

Our depth-depth-curve differs slightly from the curve which 
was gained by the visual wiggle matching (Fig. 7). How- 
ever, despite our (still) weak algorithm used to fit the non- 
parametric adjustment function to the LOS, we obtained a 
good result of adjusted data series. If they are well adjusted, 
the correlation coefficient between the parameters of the ad- 
justed data and the reference data should not vary so much. 
The correlation coefficients between the reference and ad- 
justed data series is about 0.70 - 0.80, where the correlation 
coefficients of the interactive rescaled data varies from 0.71 - 
0.87 (Tab. 1). The measure of the correlation coefficients 
emphasizes more variation for the wiggle matching than for 
the CRP re scaling. 



Depth in Core PS 21 78-3 GPC [cm] 




Fig. 7. Depth-depth-curves. In black the curve gained with the CRP, in red 
the manually matching result. The green curve shows the deviation between 
both results. 



7 Discussion 

Cross recurrence plots (CRP) reveal similarities in the states 
of the two systems. A similar trajectory evolution gives a di- 
agonal structure in the CRP. An additional time dilatation or 
compression of one of these similar trajectories causes a dis- 
tortion of this diagonal structure (Fig. 1). This effect is used 
to look into the synchronization between both systems. Syn- 
chronized systems have diagonal structures along and in the 
direction of the main diagonal in the CRP. Interruptions of 
these structures with gaps are possible because of variations 
in the amplitudes of both systems. However, a loss of syn- 
chronization is viewable by the distortion of this structures 
along the main diagonal (LOS). By fitting a non-parametric 
function to the LOS one allows to re-synchronization or ad- 
justment to both systems at the same time scale. Although 
this method is based on principles from deterministic dynam- 
ics, no assumptions about the underlying systems has to be 



Table 1. Correlation coefficients gi.2 between adjusted data and reference 
data and their deviation. The correlation of the interactive adjusted data 
varies more than the automatic adjusted data. The data length is = 170 
(wiggle matching) and A'^ = 250 (CRP matching). The difference between 
the both correlation coefficients and g2 is significant at a 99 % signifi- 
cance level, when the test measure z is greater than ^o.oi = 2.576. 



Parameter 


gi , wiggle matching 


g2 . CRP matching 


z 


ARM 


0.8667 


0.7846 


6.032 


MDFarm 


0.8566 


0.7902 


4.791 


K.LF 


0.7335 


0.7826 


2.661 


i^arm/i^lf 


0.8141 


0.8049 


0.614 


PJA 


0.7142 


0.6995 


0.675 


INC 


0.7627 


0.7966 


1.990 




141.4 


49.1 
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Fig. 8. ARM data after adjustment by wiggle matcliing (top) and by auto- 
matic adjustment using the LOS from Fig. 6. The bottom figure shows the 
reference data. 
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Depth in Core PS 2180-2 [cm] 

Fig. 9. The adjusted marine sediment parameters. The construction of the 
CRP was done with the normalized parameters. In this plots we show the 
parameters, which are not normalized. 



made in order for the method to work. 

The first example shows the obvious relationship between 
the LOS and the time domains of the considered time se- 
ries. The squared increasing of the frequency of the second 
harmonic function causes a parabolic LOS shape in the CRP 
(Fig. 2). Finally, with this LOS we are able to rescale the 
second function to the scale of the first harmonic function 
(Fig. 4). Some differences in the amplitude of the result are 
caused by the algorithm used in order to extract the LOS from 
the CRP. However, our concern is to focus on the distorted 
main diagonal and their relationship with the time domains. 

The second example deals with real geological data and al- 
lows a comparison with the result of the conventional method 
of visual wiggle matching. The visual comparison of the ad- 
justed data shows a good concordance with the reference and 
the wiggle matched data (Fig. 8 and 9). The depth-depth- 
function differs up to 20 centimeters from the depth-depth- 
function of the wiggle matching. The correlation coefficients 
between the CRP adjusted data and the reference data varies 
less than the correlation coefficients of the wiggle matching. 
However, the correlation coefficients for the CRP adjusted 
data are smaller than these for the wiggle matched data. Al- 
though their correlation is better, it seems that the interactive 
method does not produce a balanced adjusting, whereas the 
automatic matching looks for a better balanced adjusting. 

These both examples exhibits the ability to work with smooth 
and non-smooth data, whereby the result will be better for 
smooth data. Small fluctuations in the non-smooth data can 
be handled by the LOS searching algorithm. Therefore, smooth- 
ing strategies, like smoothing or parametrical fit of the LOS, 
are not necessary. The latter would damp one advantage of 
this method, that the LOS is yielded as a non-parametrical 
function. A future task will be the optimization of the LOS 
searching algorithm, in order to get a clear LOS even if the 
data are non-smooth. Further, the influence of dynamical 
noise to the result will be studied. Probably, this problem 
may be bypassed by a suitable LOS searching algorithm too. 

Our method has conspicuous similarities with the method 
of sequence slotting described by Thompson and Clark (1989). 
The first step in this method is the calculation of a distance 
matrix similar to our Eq. 2, which allows the use of multi- 
variate data sets. Thompson and Clark (1989) referred to the 
distance measure as dissimilarity. It is used to determine the 
alignment function in such a way that the sum of the dissimi- 
larities along a path in the distance matrix is minimized. This 
approach is based on dynamic programming methods which 
were mainly developed for speech pattern recognition in the 
70's (e.g. Sakoe and Chiba, 1978). In contrast, RPs were 
developed to visualize the phase space behaviour of dynam- 
ical systems. Therefore, a threshold was introduced to make 
recurrent states visible. The involving of a fixed amount of 
nearest neighbours in the phase space and the possibility to 
increase the embedding dimensions distinguish this approach 
from the sequence slotting method. 
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8 Conclusion 

The cross recurrence plot (CRP) can contain information about 
the synchronization of data series. This is revealed by the 
distorted main diagonal, which is called line of synchroniza- 
tion (LOS). After isolating this LOS from the CRP, one ob- 
tains a non-parametric rescahng function. With this function, 
one can synchronize the time series. The underlying more- 
dimensional phase space allows to include more than one pa- 
rameter in this synchronization method, as it usually appears 
in geological applications, e.g. core synchronization. The 
comparison of CRP adjusted geophysical core data with the 
conventionally visual matching shows an acceptable rehabil- 
ity level of the new method, which can be further improved 
by a better method for estimating the LOS. The advantage 
is the automatic, objective and multivariate adjustment. Fi- 
nally, this method of CRPs can open a wide range of applica- 
tions as scale adjustment, phase synchronization and pattern 
recognition for instance in geology, molecular biology and 
ecology. 
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Appendix: An Algoritlun to Fit the LOS 

In order to implement a recognition of the LOS we have used 
the following simple two-step algorithm. Denote all recur- 
rence points by Vi^j- (a, (3 — 1, 2, . . .) and the recurrence 
points lying on the LOS by (a,f3 = 1,2, . . .). Be- 

fore the common algorithm starts, find the recurrence point 
Vi^j^ next to the axes origin. In the first step, the next re- 
currence point ri~ j~ after a previous determined recurrence 
point r^^ . is to be determined. This is carried out by a step- 
wise increasing of a squared {w x w) sub-matrix, wherein the 
previous recurrence point is at the (1, 1) -location. The size 
w of this sub-matrix increases step-wise until it meets a new 
recurrence point or the margin of the CRP. When a next re- 
currence point Vi-j- = ri^+si,j/,+5j (Si = w or Sj = w) in 
the x-direction (t/-direction) is found, the second step looks 
if there are following recurrence points in y-direction (x- 
direction). If this is true (e. g. there are a cluster of recur- 
rence points) increase further the sub-matrix in j/-direction 
(x-direction) until a predefined size {w + dx) x [w + dy) 
(dx < dx, dy < dy) or until no new recurrence points are 
met. This further increasing of the sub-matrix is done for 
the both X- and y-direction. Using dx, dy we compute the 
next recurrence point rj^^j ^^^^ by determination of the cen- 
ter of mass of the cluster of recurrence points with ia+i = 
ia + {dx + 5i)/2 and = + [dy + Sj)/2. The latter 
avoids the fact that the algorithm is driven around widespread 
areas of recurrence points. Instead of this, the algorithm lo- 
cates the LOS within these areas. (However, the introduc- 
ing of two additional parameter dx and dy is a disadvan- 
tage which should be avoided in further improvements of 



this algorithm.) The next step is to set the recurrence point 
^ia+i ,3f)+i to a new start point and to begin with the step one 
in order to find the next recurrence point. These steps are 
repeated until the end of the RP is reached. 

We know that this algorithm is merely one of many possi- 
bilities. The following criteria should be met in order to ob- 
tain a good LOS. The amount of targeted recurrence points 
by the LOS A^i should converge to the maximum and the 
amount of gaps in the LOS A^o should converge to the min- 
imum. An analysis with various estimated LOS confirms this 
requirement. The correlation between two LOS-synchronized 
data series arises with Ni and with 1/N{) (the correlation co- 
efficient correlates most strongly with the ratio Ni/Nq). 

The algorithm for computation of the CRP and recogni- 
tion of the LOS are available as Matlab programmes on the 
WWW: http://www.agnld.uni-potsdam.de/~marwan. 
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